home *** CD-ROM | disk | FTP | other *** search
/ The Arsenal Files 4 / The Arsenal Files 4 (Arsenal Computer).ISO / ham / sattrk31.tgz / sattrack-3.1.tar / SatTrack / src / sattrack / satheli.c < prev    next >
C/C++ Source or Header  |  1995-03-16  |  14KB  |  297 lines
  1. /******************************************************************************/
  2. /*                                                                            */
  3. /*  Title       : satheli.c                                                   */
  4. /*  Author      : Manfred Bester                                              */
  5. /*  Date        : 17Oct93                                                     */
  6. /*  Last change : 15Mar95                                                     */
  7. /*                                                                            */
  8. /*  Synopsis    : Routines that calculate the positions of the Sun and the    */
  9. /*                Moon. If 'moonFlag' or 'sunTransitFlag' are set, the update */
  10. /*                timer is without effect.                                    */
  11. /*                                                                            */
  12. /*                                                                            */
  13. /*  SatTrack is Copyright (c) 1992, 1993, 1994, 1995 by Manfred Bester.       */
  14. /*  All Rights Reserved.                                                      */
  15. /*                                                                            */
  16. /*  Permission to use, copy, and distribute SatTrack and its documentation    */
  17. /*  in its entirety for educational, research and non-profit purposes,        */
  18. /*  without fee, and without a written agreement is hereby granted, provided  */
  19. /*  that the above copyright notice and the following three paragraphs appear */
  20. /*  in all copies. SatTrack may be modified for personal purposes, but        */
  21. /*  modified versions may NOT be distributed without prior consent of the     */
  22. /*  author.                                                                   */
  23. /*                                                                            */
  24. /*  Permission to incorporate this software into commercial products may be   */
  25. /*  obtained from the author, Dr. Manfred Bester, 1636 M. L. King Jr. Way,    */
  26. /*  Berkeley, CA 94709, USA. Note that distributing SatTrack 'bundled' in     */
  27. /*  with ANY product is considered to be a 'commercial purpose'.              */
  28. /*                                                                            */
  29. /*  IN NO EVENT SHALL THE AUTHOR BE LIABLE TO ANY PARTY FOR DIRECT, INDIRECT, */
  30. /*  SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES ARISING OUT OF THE USE OF   */
  31. /*  THIS SOFTWARE AND ITS DOCUMENTATION, EVEN IF THE AUTHOR HAS BEEN ADVISED  */
  32. /*  OF THE POSSIBILITY OF SUCH DAMAGE.                                        */
  33. /*                                                                            */
  34. /*  THE AUTHOR SPECIFICALLY DISCLAIMS ANY WARRANTIES, INCLUDING, BUT NOT      */
  35. /*  LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A   */
  36. /*  PARTICULAR PURPOSE. THE SOFTWARE PROVIDED HEREUNDER IS ON AN "AS IS"      */
  37. /*  BASIS, AND THE AUTHOR HAS NO OBLIGATIONS TO PROVIDE MAINTENANCE, SUPPORT, */
  38. /*  UPDATES, ENHANCEMENTS, OR MODIFICATIONS.                                  */
  39. /*                                                                            */
  40. /******************************************************************************/
  41.  
  42. #include <stdio.h>
  43. #include <math.h>
  44.  
  45. #ifndef STDLIB
  46. #include <stdlib.h>
  47. #endif
  48.  
  49. #include "satglobalsx.h"
  50. #include "sattrack.h"
  51.  
  52. /******************************************************************************/
  53. /*                                                                            */
  54. /* getSunVector: calculates the apparent geocentric position of the Sun       */
  55. /*               and the Moon using harmonic polynomials                      */
  56. /*                                                                            */
  57. /*               for reference see:                                           */
  58. /*               T. C. Van Flandern and K. F. Pulkkinen,                      */
  59. /*               Ap.J. (Suppl.) 41, 391-411 (1979)                            */
  60. /*                                                                            */
  61. /******************************************************************************/
  62.  
  63. void getSunVector(moonFlag)
  64.  
  65. int moonFlag;
  66.  
  67. {
  68.     double cheby[20], a[31];
  69.     double meanAnomMoon, meanLongMoon, longAscNodeMoon, argLatMoon;
  70.     double meanElongMoon, meanLongSun, meanAnomSun;
  71.     double meanAnomVenus, meanAnomMars, meanAnomJupiter;
  72.     double tu, Tu, rP, U, V, W, raa, deca;
  73.     int    i;
  74.  
  75.     if (!moonFlag && !sunTransitFlag && 
  76.         fabs(julianDate - lastJulianDateSun) < SUNUPDATEINT)
  77.         return;
  78.  
  79.     lastJulianDateSun = julianDate;
  80.  
  81.     tu = julianDate - JULDAT2000;                     /* t since 2000.0 [d]   */
  82.     Tu = tu / JULCENT + 1.0;                          /* t since 1900.0 [jcy] */
  83.  
  84.     meanLongMoon    = reduce(0.606434 + 0.03660110129 * tu, -ONE, ONE);
  85.     meanAnomMoon    = reduce(0.374897 + 0.03629164709 * tu, -ONE, ONE);
  86.     argLatMoon      = reduce(0.259091 + 0.03674819520 * tu, -ONE, ONE);
  87.     meanElongMoon   = reduce(0.827362 + 0.03386319198 * tu, -ONE, ONE);
  88.     longAscNodeMoon = reduce(0.347343 - 0.00014709391 * tu, -ONE, ONE);
  89.     meanLongSun     = reduce(0.779072 + 0.00273790931 * tu, -ONE, ONE);
  90.     meanAnomSun     = reduce(0.993126 + 0.00273777850 * tu, -ONE, ONE);
  91.     meanAnomVenus   = reduce(0.140023 + 0.00445036173 * tu, -ONE, ONE);
  92.     meanAnomMars    = reduce(0.053856 + 0.00145561327 * tu, -ONE, ONE);
  93.     meanAnomJupiter = reduce(0.056531 + 0.00023080893 * tu, -ONE, ONE);
  94.  
  95.     cheby[1]  = TWOPI * reduce(meanLongMoon,ZERO,ONE);
  96.     cheby[2]  = TWOPI * reduce(meanAnomMoon,ZERO,ONE);
  97.     cheby[3]  = TWOPI * reduce(argLatMoon,ZERO,ONE);
  98.     cheby[4]  = TWOPI * reduce(meanElongMoon,ZERO,ONE) ;
  99.     cheby[5]  = TWOPI * reduce(longAscNodeMoon,ZERO,ONE);
  100.     cheby[7]  = TWOPI * reduce(meanLongSun,ZERO,ONE);
  101.     cheby[8]  = TWOPI * reduce(meanAnomSun,ZERO,ONE);
  102.     cheby[13] = TWOPI * reduce(meanAnomVenus,ZERO,ONE);
  103.     cheby[16] = TWOPI * reduce(meanAnomMars,ZERO,ONE);
  104.     cheby[19] = TWOPI * reduce(meanAnomJupiter,ZERO,ONE);
  105.  
  106.     /* Sun */
  107.  
  108.     rP = 1.00014 - 0.01675 * cos(cheby[8]) - 0.00014 * cos(2.0 * cheby[8]);
  109.  
  110.     U =  1.00000 - 
  111.          0.03349 * cos(cheby[8]) - 
  112.          0.00014 * cos(2.0 * cheby[8]) + 
  113.          0.00008 * Tu * cos(cheby[8]) - 
  114.          0.00003 * sin(cheby[8] - cheby[19]);
  115.  
  116.     V =  0.39785 * sin(cheby[7]) - 
  117.          0.01000 * sin(cheby[7] - cheby[8]) +
  118.          0.00333 * sin(cheby[7] + cheby[8]) - 
  119.          0.00021 * Tu * sin(cheby[7]) +
  120.          0.00004 * sin(cheby[7] + 2.0 * cheby[8]) - 
  121.          0.00004 * cos(cheby[7]) -
  122.          0.00004 * sin(cheby[5] - cheby[7]) + 
  123.          0.00003 * Tu * sin(cheby[7] - cheby[8]);
  124.  
  125.     W = -0.04129 * sin(2.0 * cheby[7]) + 
  126.          0.03211 * sin(cheby[8]) +
  127.          0.00104 * sin(2.0 * cheby[7] - cheby[8]) -
  128.          0.00035 * sin(2.0 * cheby[7] + cheby[8]) - 
  129.          0.00010 -
  130.          0.00008 * Tu * sin(cheby[8]) - 
  131.          0.00008 * sin(cheby[5]) +
  132.          0.00007 * sin(2.0 * cheby[8]) + 
  133.          0.00005 * Tu * sin(2.0 * cheby[7]) +
  134.          0.00003 * sin(cheby[1] - cheby[7]) -
  135.          0.00002 * cos(cheby[8] - cheby[19]) +
  136.          0.00002 * sin(4.0 * cheby[8] - 8.0 * cheby[16] + 3.0 * cheby[19]) -
  137.          0.00002 * sin(cheby[8] - cheby[13]) -
  138.          0.00002 * cos(2.0 * cheby[8] - 2.0 * cheby[13]);
  139.  
  140.     raa    = W / sqrt(U - V*V);
  141.     sunRA  = cheby[7] + atan(raa  / sqrt(1.0 - raa*raa));
  142.     sunRA  = reduce(sunRA,ZERO,TWOPI);
  143.  
  144.     deca   = V / sqrt(U);
  145.     sunDec = atan(deca / sqrt(1.0 - deca*deca));
  146.     sunDec = reduce(sunDec,-PI,PI);
  147.  
  148.     rP    *= EARTHSMA;
  149.  
  150.     sunPosGl[0] = rP * cos(sunRA) * cos(sunDec);
  151.     sunPosGl[1] = rP * sin(sunRA) * cos(sunDec);
  152.     sunPosGl[2] = rP * sin(sunDec);
  153.  
  154.     sunDist     = absol(sunPosGl);
  155.  
  156.     if (!moonFlag)
  157.         return;
  158.  
  159.     /* Moon */
  160.  
  161.     rP    =  60.36298 
  162.             - 3.27746 * cos(cheby[2]) 
  163.             - 0.57994 * cos(cheby[2] - 2.0 * cheby[4]) 
  164.             - 0.46357 * cos(2.0 * cheby[4]) 
  165.             - 0.08904 * cos(2.0 * cheby[2]) 
  166.             + 0.03865 * cos(2.0 * cheby[2] - 2.0 * cheby[4]) 
  167.             - 0.03237 * cos(2.0 * cheby[4] - cheby[8]) 
  168.             - 0.02688 * cos(cheby[2] + 2.0 * cheby[4]) 
  169.             - 0.02358 * cos(cheby[2] - 2.0 * cheby[4] + cheby[8]) 
  170.             - 0.02030 * cos(cheby[2] - cheby[8]) 
  171.             + 0.01719 * cos(cheby[4]) 
  172.             + 0.01671 * cos(cheby[2] + cheby[8]) 
  173.             + 0.01247 * cos(cheby[2] - 2.0 * cheby[3]) 
  174.             + 0.00704 * cos(cheby[8]) 
  175.             + 0.00529 * cos(cheby[4] + cheby[8]) 
  176.             - 0.00524 * cos(cheby[2] - 4.0 * cheby[4]) 
  177.             + 0.00398 * cos(cheby[2] - 2.0 * cheby[4] - cheby[8]) 
  178.             - 0.00366 * cos(3.0 * cheby[2])
  179.             - 0.00295 * cos(2.0 * cheby[2] - 4.0 * cheby[4]) 
  180.             - 0.00263 * cos(cheby[4] + cheby[8]) 
  181.             + 0.00249 * cos(3.0 * cheby[2] - 2.0 * cheby[4]) 
  182.             - 0.00221 * cos(cheby[2] + 2.0 * cheby[4] - cheby[8]) 
  183.             + 0.00185 * cos(2.0 * cheby[3] - 2.0 * cheby[4]) 
  184.             - 0.00161 * cos(2.0 * cheby[4] - 2.0 * cheby[8]) 
  185.             + 0.00147 * cos(cheby[2] + 2.0 * cheby[3] - 2.0 * cheby[4]) 
  186.             - 0.00142 * cos(4.0 * cheby[4]) 
  187.             + 0.00139 * cos(2.0 * cheby[2] - 2.0 * cheby[4] + cheby[8]) 
  188.             - 0.00118 * cos(cheby[2] - 4.0 * cheby[4] + cheby[8])
  189.             - 0.00116 * cos(2.0 * cheby[2] + 2.0 * cheby[4]) 
  190.             - 0.00110 * cos(2.0 * cheby[2] - cheby[8]);
  191.  
  192.     a[1]  =  0.39558 * sin(cheby[3] + cheby[5]);
  193.     a[2]  =  0.08200 * sin(cheby[3]);
  194.     a[3]  =  0.03257 * sin(cheby[2] - cheby[3] - cheby[5]);
  195.     a[4]  =  0.01092 * sin(cheby[2] + cheby[3] + cheby[5]);
  196.     a[5]  =  0.00666 * sin(cheby[2] - cheby[3]);
  197.     a[6]  = -0.00644 * sin(cheby[2] + cheby[3] - 2.0 * cheby[4] + cheby[5]);
  198.     a[7]  = -0.00331 * sin(cheby[3] - 2.0 * cheby[4] + cheby[5]);
  199.     a[8]  = -0.00304 * sin(cheby[3] - 2.0 * cheby[4]);
  200.     a[9]  = -0.00240 * sin(cheby[2] - cheby[3] - 2.0 * cheby[4] - cheby[5]);
  201.     a[10] =  0.00226 * sin(cheby[2] + cheby[3]);
  202.     a[11] = -0.00108 * sin(cheby[2] + cheby[3] - 2.0 * cheby[4]);
  203.     a[12] = -0.00079 * sin(cheby[3] - cheby[5]);
  204.     a[13] =  0.00078 * sin(cheby[3] + 2.0 * cheby[4] + cheby[5]);
  205.     a[14] =  0.00066 * sin(cheby[3] + cheby[5] - cheby[8]);
  206.     a[15] = -0.00062 * sin(cheby[3] + cheby[5] + cheby[8]);
  207.     a[16] = -0.00050 * sin(cheby[2] - cheby[3] - 2.0 * cheby[4]);
  208.     a[17] =  0.00045 * sin(2.0 * cheby[2] + cheby[3] + cheby[5]);
  209.     a[18] = -0.00031 * sin(2.0 * cheby[2] + cheby[3] - 
  210.                            2.0 * cheby[4] + cheby[5]);
  211.     a[19] = -0.00027 * sin(cheby[2] + cheby[3] - 2.0 * cheby[4] + 
  212.                            cheby[5] + cheby[8]);
  213.     a[20] = -0.00024 * sin(cheby[3] - 2.0 * cheby[4] + cheby[5] + cheby[8])
  214.             -0.00021 * Tu * sin(cheby[3] + cheby[5]);
  215.  
  216.     V = 0.0;
  217.  
  218.     for (i = 1; i <= 20; i++)
  219.         V += a[i];
  220.  
  221.     U = 1.0 - 0.10828 * cos(cheby[2])
  222.             - 0.01880 * cos(cheby[2] - 2.0 * cheby[4])
  223.             - 0.01479 * cos(2.0 * cheby[4])
  224.             + 0.00181 * cos(2.0 * cheby[2] - 2.0 * cheby[4])
  225.             - 0.00147 * cos(2.0 * cheby[2])
  226.             - 0.00105 * cos(2.0 * cheby[4] - cheby[8])
  227.             - 0.00075 * cos(cheby[2] - 2.0 * cheby[4] + cheby[8])
  228.             - 0.00067 * cos(cheby[2] - cheby[8])
  229.             + 0.00057 * cos(cheby[4])
  230.             + 0.00055 * cos(cheby[2] + cheby[8])
  231.             - 0.00046 * cos(cheby[2] + 2.0 * cheby[4])
  232.             + 0.00041 * cos(cheby[2] - 2.0 * cheby[3])
  233.             + 0.00024 * cos(cheby[8]);
  234.  
  235.     a[1]  =  0.10478 * sin(cheby[2]);
  236.     a[2]  = -0.04105 * sin(2.0 * cheby[3] + 2.0 * cheby[5]);
  237.     a[3]  = -0.02130 * sin(cheby[2] - 2.0 * cheby[4]);
  238.     a[4]  = -0.01779 * sin(2.0 * cheby[3] + cheby[5]);
  239.     a[5]  =  0.01774 * sin(cheby[5]);
  240.     a[6]  =  0.00987 * sin(2.0 * cheby[4]);
  241.     a[7]  = -0.00338 * sin(cheby[2] - 2.0 * cheby[3] - 2.0 * cheby[5]);
  242.     a[8]  = -0.00309 * sin(cheby[8]);
  243.     a[9]  = -0.00190 * sin(2.0 * cheby[3]);
  244.     a[10] = -0.00144 * sin(cheby[2] + cheby[5]);
  245.     a[11] = -0.00144 * sin(cheby[2] - 2.0 * cheby[3] - cheby[5]);
  246.     a[12] = -0.00113 * sin(cheby[2] + 2.0 * cheby[3] + 2.0 * cheby[5]);
  247.     a[13] = -0.00094 * sin(cheby[2] - 2.0 * cheby[4] + cheby[8]);
  248.     a[14] = -0.00092 * sin(2.0 * cheby[2] - 2.0 * cheby[4]);
  249.     a[15] =  0.00071 * sin(2.0 * cheby[4] - cheby[8]);
  250.     a[16] =  0.00070 * sin(2.0 * cheby[2]);
  251.     a[17] =  0.00067 * sin(cheby[2]+2.0 * cheby[3] -
  252.                            2.0 * cheby[4]+2.0 * cheby[5]);
  253.     a[18] =  0.00066 * sin(2.0 * cheby[3] - 2.0 * cheby[4] + cheby[5]);
  254.     a[19] = -0.00066 * sin(2.0 * cheby[4] + cheby[5]);
  255.     a[20] =  0.00061 * sin(cheby[2] - cheby[8]);
  256.     a[21] = -0.00058 * sin(cheby[4]);
  257.     a[22] = -0.00049 * sin(cheby[2] + 2.0 * cheby[3] + cheby[5]);
  258.     a[23] = -0.00049 * sin(cheby[2] - cheby[5]);
  259.     a[24] = -0.00042 * sin(cheby[2] + cheby[8]);
  260.     a[25] =  0.00034 * sin(2.0 * cheby[3] - 2.0 * cheby[4] + 2.0 * cheby[5]);
  261.     a[26] = -0.00026 * sin(2.0 * cheby[3] - 2.0 * cheby[4]);
  262.     a[27] =  0.00025 * sin(cheby[2] - 2.0 * (cheby[3] + cheby[4] + cheby[5]));
  263.     a[28] =  0.00024 * sin(cheby[2] - 2.0 * cheby[3]);
  264.     a[29] =  0.00023 * sin(cheby[2] + 2.0 * cheby[3] -
  265.                                       2.0 * cheby[4] + cheby[5]);
  266.     a[30] =  0.00023 * sin(cheby[2] - 2.0 * cheby[4] - cheby[5]);
  267.  
  268.     W = 0.0;
  269.  
  270.     for (i = 1; i <= 30; i++)
  271.         W += a[i];
  272.  
  273.     raa     = W / sqrt(U - V*V);
  274.     moonRA  = cheby[1] + atan(raa  / sqrt(1.0 - raa*raa));
  275.     moonRA  = reduce(moonRA,ZERO,TWOPI);
  276.  
  277.     deca    = V / sqrt(U);
  278.     moonDec = atan(deca / sqrt(1.0 - deca*deca));
  279.     moonDec = reduce(moonDec,-PI,PI);
  280.  
  281.     rP     *= EARTHRADIUS;
  282.  
  283.     moonPosGl[0] = rP * cos(moonRA) * cos(moonDec);
  284.     moonPosGl[1] = rP * sin(moonRA) * cos(moonDec);
  285.     moonPosGl[2] = rP * sin(moonDec);
  286.  
  287.     moonDist     = absol(moonPosGl);
  288.  
  289.     return;
  290. }
  291.  
  292. /******************************************************************************/
  293. /*                                                                            */
  294. /* End of function block satheli.c                                            */
  295. /*                                                                            */
  296. /******************************************************************************/
  297.